%% ============ formal codes for fracture generatio from seismicity events==========
function [frac_info, loss_goal,fit] = fracture_generation(data,seis)
    strike_min=seis.s_min;
    strike_max=seis.s_max;
    num_frac=seis.num_frac;
    threshold=seis.buffer;
    loss_goal=seis.loss_;
    
    
    n_angle=size(strike_min,2);
    
    fixed_=seis.fix_fracs;
    fixed=fixed_;
    nf=size(fixed_,1);
    idx=randperm(nf);
    for i =1:nf
        fixed(i,:)=fixed_(idx(i),:);
    end

    w=seis.weight;

    %=============judge if any fixed fractures =========================
    if  isempty(fixed)
        stochast=0;
        fixed=[];
    else
        stochast=1;
    end

    for i = 1:seis.iters
        lines = 0  ;
        f_info = [];
        data_update=data;
        fitted1= [];
        goal_dist = 0;
        
        while length(data_update) > length(data) * 0.05

            len_data = length(data_update);
            rcenter = data_update(randi(len_data),:)  ;
            ran_number = randi(n_angle) ;      
            
            %strike = round(unifrnd(strike_min(ran_number), strike_max(ran_number)));
            strike = unifrnd(strike_min(ran_number), strike_max(ran_number));
            
            if lines+1<=size(fixed,1) && stochast>0.5
                rcenter = fixed(lines+1,1:2) ;
                strike = fixed(lines+1,3);
            end

            rcenter=[rcenter,strike];

            fitted = [];
            positions=[];

            for dat=1:len_data
                [dist, propit] = cal_distance(rcenter, data_update(dat,:));
                goal_dist = goal_dist+dist;
                if dist <= threshold 
                    fitted    = [fitted; data_update(dat,:)];
                    fitted1   = [fitted1; data_update(dat,:)];
                    positions = [positions; propit];
                end
            end
            data_update=setdiff(data_update,fitted,'rows');
            if strike<90
                startX = min(positions(:,1));
                endX   = max(positions(:,1));
                startY = min(positions(:,2));
                endY   = max(positions(:,2)); 
            else
                startX = min(positions(:,1));
                endX   = max(positions(:,1));
                startY = max(positions(:,2));
                endY   = min(positions(:,2)); 
            end
            frac_length=sqrt((startX-endX)^2+(startY-endY)^2) ;
            
            if frac_length>2.0
                frac_i=[startX,endX,startY,endY, frac_length, strike, seis.aperture];
                f_info = [f_info; frac_i];
                lines = lines+ 1;
            end

            if lines >= num_frac
                break
            end
        end
        
        fit_ratio=1-length(data_update)/length(data);
        fit_target=goal_dist/w/fit_ratio;
        
        loss=fit_target;
    
        if loss<loss_goal
            loss_goal     = loss;
            frac_info     = f_info;
            fit=fitted1;
            % N_fracture = lines;
        end 
      
    end
    

end
